#include "./lj-cutoff.h"
#include <math.h>

void mds_lj_cutoff_init( mds_lj_cutoff_data *lj_cutoff, Real epsilon, Real sigma, Real r_c )
{
    assert(lj_cutoff);
	memset( lj_cutoff, 0, sizeof(mds_lj_cutoff_data) );
	fprintf( stderr, "  [Initializing LJ cutoff: epsilon=%g sigma=%g r_c=%g]\n", epsilon, sigma, r_c);
    lj_cutoff->epsilon   = epsilon;                /*!< in kJ/mol */
    lj_cutoff->epsilon4  = MDS_F(4.0)  * epsilon;  /*!< 4 epsilon */
    lj_cutoff->epsilon24 = MDS_F(24.0) * epsilon;  /*!< 24 espilon */
    lj_cutoff->sigma     = sigma;                  /*!< in nm     */
    lj_cutoff->r_c       = r_c;
    lj_cutoff->r_csq     = r_c * r_c;
    lj_cutoff->sigma_sq  = sigma * sigma;   /*!< sigma^2   */ 
}

Real mds_lj_cutoff( const mds_vector *r, mds_vector *f, void *args )
{
    assert( r );
    assert( f );
    assert( args );
	{
		const mds_lj_cutoff_data *lj_cutoff = args;
        const Real   r2       				= mds_vector_normsq(r);
        Real		 r1						= (Real)sqrt( r2 );
        const Real   r_c2    				= lj_cutoff->r_csq;
        const Real   q2       				= lj_cutoff->sigma_sq / r2;
		const Real	 q2_c     				= lj_cutoff->sigma_sq / r_c2;
        const Real   q6       				= q2 * q2 * q2;
		const Real	 q6_c     				= q2_c * q2_c * q2_c;
        const Real   q12      				= q6 * q6;
		const Real   q12_c    				= q6_c * q6_c;
        const Real   q12mq6   				= q12 - q6;
		const Real   q12mq6_c 				= q12_c - q6_c;
		
		if( r2 <= lj_cutoff->r_csq )
		{
			/** compute force **/
			mds_vector_mulset(f, lj_cutoff->epsilon24 * (( q12 + q12mq6 ) / r2 - ((1/r1) * (1/lj_cutoff->r_c) * (q12_c + q12mq6_c))), r);
			
			/** compute potential */
			return (lj_cutoff->epsilon4 * ( q12mq6 - q12mq6_c ) + ((r1-lj_cutoff->r_c)/lj_cutoff->r_c) * lj_cutoff->epsilon24 * (q12_c + q12mq6_c)); 
		}
		else                                      
		{
			mds_vector_ldz(f); /*return \vec{f} = 0*/ 
			return 0;
		}
	}
}

#include "./lj.h"
void mds_ljc_init( mds_ljc_data *ljc, Real  epsilon, Real sigma, Real r_c )
{
	assert(ljc);
	const mds_vector   R = { r_c, 0, 0 };
	mds_vector         F = { 0, 0, 0 };
	ljc->sigma_sq        = sigma * sigma;
	ljc->epsilon4        = MDS_F(4.0)  * epsilon;
    ljc->epsilon24       = MDS_F(24.0) * epsilon;
	ljc->r_c             = r_c;
	ljc->r_c_sq          = r_c * r_c;
	ljc->V_c             = mds_lj( &R, &F, ljc );
	ljc->F_c             = F.x;
	
}

Real mds_ljc( const mds_vector *r, mds_vector *f, void *args )
{
	assert( r );
	assert( f );
	assert( args );
	
	{
		const mds_ljc_data *ljc = args;
		const Real          r2  = mds_vector_normsq(r);
		if( r2 < ljc->r_c_sq )
		{
			const Real         q2     = ljc->sigma_sq / r2; 
			const Real         q6     =  q2 * q2 * q2;
			const Real         q12    =  q6 * q6;
			const Real         r1     = (Real)sqrt(r2);
			
			/** compute force **/
			mds_vector_mulset(f, ljc->epsilon24 * ( (q12+q12) - q6 ) / r2 - ljc->F_c/r1, r);
			
			/** return potential **/
			return ( ( ljc->epsilon4 * ( q12 -  q6 ) ) - ljc->V_c ) + (r1-ljc->r_c) * ljc->F_c;
		}
		else 
		{
			mds_vector_ldz( f );
			return 0;
		}
		
	}
	
}
